InĀ [1]:
set.seed(999)
options(scipen = 9)
options(warn = -1) 
source("./environment/libraries.R")
knitr::opts_chunk$set(fig.height = 12, fig.width = 18, fig.dpi = 300)
knitr::opts_chunk$set(warning = FALSE)
All packages loaded successfully
InĀ [2]:
name <- "Kenya_E1"
dataset <- read.csv(paste0("./test/", name, "_processed.csv"))
dataset_c <- data.frame(dataset[7:ncol(dataset)], row.names = dataset$Serial)
proj_coord <- data.frame("Easting" = dataset$Easting, "Northing" = dataset$Northing)
xy_coord <- data.frame("Longitude" = dataset$Longitude, "Latitude" = dataset$Latitude)

dataset_c_closed <- cbind(dataset_c, "Res" = 100 - rowSums(dataset_c)) 
head(dataset_c_closed)

dataset_spdf <- cbind(proj_coord, dataset_c_closed)
rownames(dataset_spdf) <- rownames(dataset_c)

structures_geojson_path <- file.path("./data", paste0(name, "_topo_lines.geojson"))
structures <- st_read(structures_geojson_path, quiet = TRUE)
structures_utm <- st_transform(structures, crs = 32637) # Transform to UTM Zone 35S (WGS84)


output_dirs <- "./test/ck/"
if (!dir.exists(output_dirs)) {
  dir.create(output_dirs, recursive = TRUE)
}
A data.frame: 6 Ɨ 17
PCaTiMnFeNiZnRbSrYZrBaMgAlSiKRes
<dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl>
10.27313.47950.33090.06032.72480.00210.00330.00350.08080.00170.0123000000.09331.2275275.45467821.960071.97141662.32071
20.38353.95870.28280.05662.40540.00170.00310.00370.08390.00160.0041000000.09551.2341004.31907318.230192.22713966.70890
30.41584.07200.22650.06172.16340.00110.00330.00340.09210.00120.0045476330.07551.3168314.05277517.194472.48983467.82554
40.37043.95890.29020.06012.63740.00220.00350.00300.08370.00210.0102000000.08801.3842474.66847818.721382.42592465.29028
50.29813.39790.29840.07072.70030.00260.00300.00350.07970.00170.0152000000.10051.2800484.96270820.899372.48651463.39976
60.43734.00800.18080.05371.92700.00100.00290.00300.08600.00120.0035276360.05831.3045303.23393613.664892.97314372.06078
InĀ [3]:
source("./utils/functions/create_quick_map.R")
options(repr.plot.width = 10, repr.plot.height = 10)

dbscan_result <- dbscan(dataset_spdf[, c("Easting", "Northing")], eps = 10, minPts = 5)
Area <- factor(dbscan_result$cluster)
dataset$Area <- Area
create_quick_map(dataset, structures, group_data = "Area")

dataset_spdf$Area <- Area
dataset_spdf <- dataset_spdf[dataset_spdf$Area == 1, ]
No description has been provided for this image
InĀ [4]:
create_spatial_objects <- function(data_spdf, transformation = "none") {
  coords <- data_spdf[, c("Easting", "Northing")]
  crs_string <- "+proj=utm +zone=37 +north +datum=WGS84"
  
  if (transformation == "ilr") {
    spatial_data <- as.data.frame(ilr(data_spdf[, -c(1, 2, ncol(data_spdf))]))
  } else if (transformation == "clr") {
    spatial_data <- as.data.frame(clr(data_spdf[, -c(1, 2, ncol(data_spdf))]))
  } else {
    spatial_data <- data_spdf[, -c(1, 2, ncol(data_spdf))]
  }
  
  return(SpatialPointsDataFrame(coords = coords, data = spatial_data, 
                               proj4string = CRS(crs_string)))
}

spdf_comp <- create_spatial_objects(dataset_spdf, "none")
spdf_ilr <- create_spatial_objects(dataset_spdf, "ilr")
spdf_clr <- create_spatial_objects(dataset_spdf, "clr")

source("./utils/functions/lag_distance_from_spdf.R")
source("./utils/functions/site_diagonal_from_spdf.R")
lag_dist <- lag_distance_from_spdf(spdf_ilr)
site_diag <- site_diagonal_from_spdf(spdf_ilr)
InĀ [5]:
source("./utils/functions/geostatistical_workflow.R")
InĀ [6]:
results_ilr_omni <- geostatistical_workflow(spdf_ilr, "ILR_Omnidirectional", lag_dist, site_diag, directional = FALSE)
Selected models:
         psill    range kappa   n
Nug 0.00105033 0.000000     0 136
Gau 0.02873073 3.639113     0  78
Exp 0.01735112 3.593177     0  58
Linear Model of Coregionalization found. Good.
[using ordinary cokriging]
Saved 17 maps to: ./test/ck/ilr-omnidirectional/maps/ 
                V1             V2             V3             V4             V5
ME  -0.00002187111  0.02407150016  0.00281125028 -0.00226650633 -0.00243363142
MSE     0.04929753     0.09519844     0.02218845     0.02013800     0.03609853
                V6             V7             V8             V9            V10
ME   0.00049795368  0.00378105592 -0.00102092507  0.00646904283  0.00208785007
MSE     0.01606958     0.01892075     0.01062964     0.02852003     0.20423875
               V11            V12            V13            V14            V15
ME   0.00311281970 -0.00211847017  0.00763336337  0.00752284718 -0.00555857358
MSE     0.02407379     0.01572140     0.02125674     0.02344569     0.08470890
               V16
ME  -0.00015847540
MSE     0.01371514
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
InĀ [7]:
results_ilr_dir <- geostatistical_workflow(spdf_ilr, "ILR_Directional", lag_dist, site_diag, directional = FALSE)
Selected models:
         psill    range kappa   n
Nug 0.00105033 0.000000     0 136
Gau 0.02873073 3.639113     0  78
Exp 0.01735112 3.593177     0  58
Linear Model of Coregionalization found. Good.
[using ordinary cokriging]
Saved 17 maps to: ./test/ck/ilr-directional/maps/ 
               V1            V2            V3            V4            V5
ME  -0.0143490895 -0.0104323151 -0.0081959822 -0.0083762146 -0.0182540391
MSE   0.048068220   0.087240677   0.020802082   0.019074672   0.032819468
               V6            V7            V8            V9           V10
ME   0.0020152234 -0.0007858862  0.0008997367 -0.0148447382 -0.0304175161
MSE   0.014418483   0.018582763   0.009193277   0.026394854   0.213054365
              V11           V12           V13           V14           V15
ME  -0.0086429833  0.0055115959 -0.0094178968 -0.0103814351  0.0244454560
MSE   0.021329791   0.014623150   0.023473660   0.026800473   0.082047630
              V16
ME   0.0049010919
MSE   0.011445260
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
InĀ [8]:
par(bg = "white")
options(repr.plot.width = 15, repr.plot.height = 12)

# Compute MAF analysis
g_temp <- create_gstat_from_spdf(spdf_ilr, method = "ordinary")
v_omni <- variogram(g_temp, width = lag_dist / 2, cutoff = site_diag / 3, cross = TRUE)
maf_result <- Maf(spdf_ilr@data, vg = v_omni)

source("./utils/functions/quick_pca_screeplot.R")
quick_pca_screeplot(maf_result)
maf_dataset <- maf_result$scores[, 1:9]
spdf_maf <- SpatialPointsDataFrame(coords = dataset_spdf[, c("Easting", "Northing")], 
                                  data = as.data.frame(maf_dataset),
                                  proj4string = CRS("+proj=utm +zone=37 +north +datum=WGS84"))

results_maf_omni <- geostatistical_workflow(spdf_maf, "MAF_Omnidirectional", lag_dist, site_diag,
                                           directional = FALSE, maf_result = maf_result, 
                                           resolution = 2)
Selected models:
        psill     range kappa  n
Nug 0.0000000 0.0000000     0 45
Exp 1.0927552 0.9507579     0 34
Gau 0.1560802 9.9853154     0 11
Linear Model of Coregionalization found. Good.
[using ordinary cokriging]
No description has been provided for this image
Saved 17 maps to: ./test/ck/maf-omnidirectional/maps/ 
            maf1         maf2         maf3         maf4         maf5
ME   0.003786155 -0.023257962  0.014695922 -0.040568633  0.000701461
MSE    1.0337736    0.9979298    1.0528312    1.3033154    1.3210178
            maf6         maf7         maf8         maf9
ME  -0.005600492  0.013122820 -0.038253927 -0.031507740
MSE    1.2730757    1.1258029    0.8005447    1.2387432
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
InĀ [9]:
g_temp2 <- create_gstat_from_spdf(spdf_ilr, method = "ordinary")
v_dir <- variogram(g_temp2, width = lag_dist / 2, cutoff = site_diag / 3,
                  alpha = c(0, 45, 90, 135), tol.hor = 22.5, cross = FALSE)

maf_result2 <- Maf(spdf_ilr@data, vg = v_dir)
quick_pca_screeplot(maf_result2)
maf_dataset2 <- maf_result2$scores[, 1:10]
spdf_maf2 <- SpatialPointsDataFrame(coords = dataset_spdf[, c("Easting", "Northing")], 
                                   data = as.data.frame(maf_dataset2),
                                   proj4string = CRS("+proj=utm +zone=37 +north +datum=WGS84"))

results_maf_dir <- geostatistical_workflow(spdf_maf2, "MAF_Directional", lag_dist, site_diag,
                                          directional = TRUE, maf_result = maf_result2, 
                                          resolution = 2)
Directional approach not fully implemented yet - using omnidirectional as fallback
Selected models:
        psill     range kappa  n
Nug 0.9821561  0.000000     0 55
Exp 0.1721345  2.481358     0 35
Gau 0.2113754 13.329098     0 12
Sph 0.4112454  2.030358     0  8
Linear Model of Coregionalization found. Good.
[using ordinary cokriging]
No description has been provided for this image
Saved 17 maps to: ./test/ck/maf-directional/maps/ 
            maf1         maf2         maf3         maf4         maf5
ME  -0.006610265  0.012043116 -0.034670502 -0.006719259 -0.018199105
MSE    1.0824425    1.1081471    1.0520101    0.9895489    1.1084641
            maf6         maf7         maf8         maf9        maf10
ME   0.038734058 -0.048576979 -0.016445294  0.082279189 -0.047689524
MSE    1.0972327    1.1123020    0.7993823    1.1104866    0.9508502
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
InĀ [10]:
source("./utils/functions/create_ck_maplist_from_list.R")

ck_files <- list.files("./test/ck/", full.names = TRUE, pattern = paste0(name, ".*\\.rds$"), recursive = TRUE)
ck_list <- lapply(ck_files, readRDS)

if (length(ck_list) > 0) {
  cat("Available cokriging results:", length(ck_list), "\n")
  cat("Files:", basename(ck_files), "\n")
  
  if (length(ck_list) > 1) {
    cat("Multiple methods completed. Individual results saved in method subdirectories\n")
  }
} else {
  cat("No cokriging results found\n")
}
cat("Geostatistical workflow completed successfully!\n")
Available cokriging results: 4 
Files: Kenya_E1_CK1.rds Kenya_E1_CK1.rds Kenya_E1_CK1.rds Kenya_E1_CK1.rds 
Multiple methods completed. Individual results saved in method subdirectories
Geostatistical workflow completed successfully!